Monte Carlo simulation of boson lattices 
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Boson lattices are theoretically well described by the Hubbard model. The basic model and its 
variants can be effectively simulated using Monte Carlo techniques. We describe two newly developed 
approaches, the Stochastic Series Expansion (SSE) with directed loop updates and continuous-time 
Diffusion Monte Carlo (CTDMC). SSE is a formulation of the finite temperature partition function 
as a stochastic sampling over product terms. Directed loops is a general framework to implement 
this stochastic sampling in a non-local fashion while maintaining detailed balance. CTDMC is 
well suited to finding exact ground-state properties, applicable to any lattice model not suffering 
from the sign problem; for a lattice model the evolution of the wave function can be performed 
in continuous time without any time discretization error. Both the directed loop algorithm and 
the CTDMC are important recent advances in development of computational methods. Here we 
present results for a Hubbard model for anti-ferromagnetic spin-1 bosons in one dimensions, and 
show evidence for a dimerized ground state in the lowest Mott lobe. 
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I. INTRODUCTION 

An optical lattice can be made by applying orthogonal laser beams to an ultracold gas of atoms. As a result 
*^Rb or ^"^Na atoms can be trapped to form a perfect lattice. These atoms have a hyperfine spin 1, with either 
a ferromagnetic (^^Rb) or antiferromagnetic (^^Na) interaction. Unpolarized ^"^Na atoms have spin-correlated Mott 
insulating states^ , and it has been suggested^ that the ground state is a dimer phase in one, two and three dimensions. 
An effective Hamiltonian of spin-1 bosons in an optical lattice has the Bose-Hubbard form, supplemented with a term 
that describes the spin interaction of atoms on the same lattice site^, 

H = -tY, («L%a + cc) - A^E"' + I - 1) + Y 51 (^^ - 2n,) , (1.1) 
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where the spin components a = x,y,z form a basis where the spin operator at site i is written in terms of boson 
operators aj, as S" — — ze^'^'^ajjO-y. (e"^"' is the totally antisymmetric Levi-Civita tensor). In this basis the antiferro- 
magnetic spin interaction has no sign problem. The value of U2/U can be determined from experimental scattering 
lengths. 

SSEi^ is an exact scheme, where the quantum partition function is expanded as a power series in a given basis {js}}, 

s n— 

where (3 — \/{kBT) is the inverse temperature. Typically the low-temperature phases we are simulating have at 
most 1-3 atoms per site, so even with the spin degrees of freedom the required number of single-site states is rather 
small. The Hamiltonian (|1.1|) couples at most two neighboring sites, so SSE splits the partition function to a sum of 
products of bond and site operations. 

To generate the terms in SSE we employ directed loop updates^. The basic idea is to pick a few well chosen 
elementary update operations, that change the bond or site operators to each other. Each update affects only a 
single site (and operator) at a time and is applied in a loop among the operators in the product that is the current 
term of the SSE sum. The outcome is a new product of operators, and, thanks to the looping, this new product is 
also a term in the partition function. In addition one controls the number of operators in the product by adding or 
removing diagonal operators; they don't change the state so their addition or removal won't disrupt the continuity of 
the sequence. The rules for how the loop travels follow from the detailed balance condition, which is used to ascertain 
that the operators appear in the products with proper weights. 

If the Hilbert space is finite one can formulate Diffusion Monte Carlo (DMC) in continuous timeSii. DMC is a 
stochastic simulation of the imaginary-time evolution operator e"^"^. For an infinitesimal time step dr this operator 
describes one out of three possible actions: 1) No action, the current state is unchanged, 2) Transition to a new state. 
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and 3) The weight of the current state changes. Action 3) is needed due to the in general non-Markovian nature 
of the evolution operator. The probabilities for the different actions p{i), i = 1,2,3 can be read directly off the 
Hamiltonian. The probability p(2) is determined by off-diagonal elements of the Hamiltonian and is therefore of order 
dr. Also p{3) is of order dr because it describes the deviation from Markovian evolution, thus only p{l) is of order 
unity. It follows that the DMC can be treated as a continuous-time simulation of a radioactive (multi-channel) decay 
problem, where decay-times are generated according to the exponential distribution: rdccay = — ln{r)dT/{l ~ p{l)) 
where r is a random number uniformly distributed between and 1. Having generated a decay-time the system is 
moved directly to the time of decay and the appropriate decay process is chosen dependent on the ratio p{2)/p{3). 
CTDMC can be combined with other standard improvements of DMC such as importance sampling, reweighting and 
forward- walking^. 



II. RESULTS 




FIG. 1: The phase boundaries between Mott insulating and superfluid regions of a ID chain for U2/U = 0.2 (left panel) and 
U2/U = 0.4 (right panel) at zero temperature. The results were obtained using CTDMC to compute ground state energies 
El{N) for different particle numbers N and system sizes L. The upper(lower) boundary /ip+(/ip_) of the Mott lobe with p 
bosons per site was obtained from the finite size values for fip± = ± {El{pL ± 1) — El{pL)) extrapolated to L ^ cxa using 
L = 8, 16, 20 and 24. 

To study the structure of the Mott insulating phase with one atom per site we have measured the bond hopping or 
dimer susceptibility 

X{q) = ^((^1^-^) ~ WiN-,)) , (2.1) 

where Ng is the Fourier transform of iVj, the number of hopping operators on bond i. Because of the on-site spin 
scattering, the hopping operators tend to couple pairwise adjacent sites, indicated by a peak at wave number q — n. 
For a dimer state in an infinite system this peak height should diverge, and we demonstrate this in Fig. |21 For 
larger than 32 sites the computation of the susceptibility becomes exceedingly slow, as autocorrelation times increase 
rapidly. This is due to our non-optimal treatment of the terms of type 

'^3c 

which changes two spin indices on 

the same site simultaneously. In order to build such terms out of loop changes where only one spin index changes 
we have included intermediate, auxiliary terms aj,ay in the Hamiltonian. These auxiliary terms are generally present 
when the loop is being built, but we do not allow a loop to close if any of them are present. 

The left panel of Fig. |3| show how the dimer state is suppressed as one moves from the insulating to the super- 
fluid state. The hopping parameter t is kept fixed and we move out of the insulating phase by increasing /i. The 
corresponding density, as it increases from unity, is shown in the right panel of Fig. |3| 
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FIG. 2: Finite size scaling of the dimer susceptibility X{q) at q = tt in the Mott insulating phase of a ID chain at /3 = 150 
for two values of U2/U indicated in the figure. We show the results for L = 4, 8, 16 and 32 sites, those for U2/U = 0.4 were 
computed at t = 0.1, /i — 0.1, and those for U2/U — 0.2 a.t t — 0.15, fi = 0.25. 




FIG. 3: The left panel shows the dimer susceptibility X{q) at the upper edge of the first Mott lobe (see the left panel of Fig. 1) 
of a ID chain with 32 sites. The susceptibility is plotted as a function of the chemical potential n/U and the wave number q. 
The data was computed at t = 0.15, l3 = 150 and U2/U = 0.2. The statistical error is less than 0.2 in the vertical scale. The 
right panel shows the corresponding density. 
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